	SUBROUTINE FMATD(X,Y,XM,YM,DB,FI,DFI,alfa,KODE,INC,NG,NNP,NP,
     &			NFIX,NN,N_eq)

	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	include 'param_dim.inc'

	DIMENSION X(NX),Y(NX),XM(NX),YM(NX),FI(NX),DFI(NX),
     *          DH(NX,NX),DG(NX,NX),DB(N_eq,NNP*2),NP(NX)
	DIMENSION INC(NX,2),KODE(NX)
	COMMON /incdata/ C0,CINC,NNINC,NUMPI,NUMPCI,NINCL,VARIA

	DIMENSION NNINC(ninc_max),CINC(ninc_max)
	integer NUMPI(ninc_max),NUMPCI(ninc_max)
	integer VARIA(ninc_max)
	DIMENSION PGLIM(nparg_max,3)

	DO 10 IP=1,NNP !loop dos nohs parametros


		NEPP= NP(IP) !numero do elemento posterior ao noh parametro
		NEAP= NEPP -1 !numero do elemento anterior ao noh parametro


		if(INC(NEAP,2).ne.NP(IP))then  ! procurar NEAP se o noh parametro for o 1o de cada spline
			do j=NFIX+1,NN
				if(INC(j,2).eq.NP(IP))then
					NEAP = j
				endif
			enddo
		endif



		do j=1,N_eq  !inicialização das matrizes DH e DG
			do i=1,N_eq
				DG(i,j)=0.
				DH(i,j)=0.
			enddo
		enddo



		DO 20 ID=1,2   ! derivada nas duas direcoes

			IC= (IP-1)*2 + ID !coluna da matriz DB referente ao noh parmetro IP, direcao ID

			DO 30 I=1,NN  ! loop sobre os pontos de colocacao (fixos e moveis)
				DO 40 J=1,NN ! loop sobre os elementos (fixos e moveis)
					INFA= 0 ! ponto de colocacao em elemento alterado INFA = 1
					INCA= 0 ! integracao em elemento alterado INCA = 1
					IF (I.EQ.NEPP) THEN  ! ponto de colocacao em elemento alterado pelo noh parametro

						INFA= 1
						IF (I.EQ.J) THEN
							CALL INLOD(X(inc(J,1)),Y(inc(J,1)),
     &								   X(inc(J,2)),Y(inc(J,2)),
     *								   DG(I,J),NEAP,ID,J,alfa)
							DH(I,J)= 0.
						ELSE IF ( J.EQ.NEAP) THEN
							INCA= 1
							CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &									       X(inc(J,2)),Y(inc(J,2)),
     *							DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)
						ELSE
							CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &									       X(inc(J,2)),Y(inc(J,2)),
     *  						DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)
                     				END IF

					ELSE IF (I.EQ.NEAP) THEN !ponto de colocacao em elemento alterado pelo noh parametro

                    				INFA= 1
						IF (I.EQ.J) THEN
							CALL INLOD(X(inc(J,1)),Y(inc(J,1)),
     &								   X(inc(J,2)),Y(inc(J,2)),
     *								DG(I,J),NEAP,ID,J,alfa)
							DH(I,J)= 0.
						ELSE IF ( J.EQ.NEPP) THEN
							INCA= 1
							CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &									       X(inc(J,2)),Y(inc(J,2)),
     *							DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)

						ELSE
							CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &									       X(inc(J,2)),Y(inc(J,2)),
     *							DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)
						END IF

					ELSE IF ((J.EQ.NEAP).OR.(J.EQ.NEPP)) THEN !ponto de colocacao em elemento NAO alterado e integracao sobre elementos alterados
						INCA= 1
						CALL INTED(XM(I),YM(I),X(inc(J,1)),Y(inc(J,1)),
     &								       X(inc(J,2)),Y(inc(J,2)),
     *  					DH(I,J),DG(I,J),NG,INCA,INFA,NEAP,ID,J,alfa)

					END IF

 40            			CONTINUE
 30         		CONTINUE

! compatibilidades nas inclusoes (analogo a montagem das matrizes H e G)

			nna = NN ! inicializa o numero de nos acumulado.
			nna0= NFIX

			do iinc= 1,nincl
				nni= nninc(iinc)  !Numero de nos (elementos) da inclusao

				do j=1,nni
					do i = 1,NN ! linha de DH ate NN
						DH(i,nna+j)= -DG(i,nna0+j)!copia as colunas de DG referentes ao fluxo nos elementos das inclusoes
					end do

					do i = 1,nni ! linhas de DH de NN ate N_eq
						DH(nna+i,nna+j)= c0*DG(nna0+i,nna0+j)/cinc(iinc) ! coloca em DH as linhas de DG afetadas da relacao entre condutividades
						if (i.eq.j) then
							DH(nna+i,nna0+j)=  DH(nna0+i,nna0+j) !DH eh sempre zero na diagonal
						else
							DH(nna+i,nna0+j)= -DH(nna0+i,nna0+j) !copia parte de DH com sinal invertido
						endif
					end do
				end do

				nna = nna  + nni ! atualiza o numero de nohs acumulados
				nna0= nna0 + nni
			enddo !fim do loop da relacao de compatibilidades


			DO 60 J=1,NFIX !Rearranjo das matrizes DH(elemntos referentes as incog ja calc DFI) e DG(valores prescritos FI)
				IF(KODE(J).EQ.0) THEN
					DO 70 I=1,NN !ate NN pq de NN ate N_eq 
						CH= DH(I,J)
						DH(I,J)= -DG(I,J)
						DG(I,J)= -CH
 70               			CONTINUE
				END IF
 60         		CONTINUE


			do i=1,N_eq
				DB(I,IC)= 0.
			enddo


			DO J=1,NFIX ! FI possui apenas NFIX linhas
				DO I=1,N_eq ! Montagem da coluna IC
					DB(I,IC)= DB(I,IC) + DG(I,J)*FI(J)
				enddo
			ENDDO
			DO J=1,N_eq ! DFI possui N_eq linhas
				DO I=1,N_eq ! Montagem da coluna IC
					DB(I,IC)= DB(I,IC) - DH(I,J)*DFI(J)
				enddo				
			ENDDO


 20     	CONTINUE !fim do loop das direcoes
 10   	CONTINUE !fim do loop dos nohs parametros 



      RETURN
      END

